The growth of galaxies in cosmological simulations of structure formation 



Chigurupati Murali, Neal Katz 
Astronomy Department, University of Massachusetts, Amherst, MA 01003 

Lars Hernquist 

■ Department of Astronomy, Harvard University, Cambridge, MA 02138 

o ■ 

psj ' David H. Weinberg 

c ■ 

^ ' Astronomy Department, Ohio State University, Columbus, OH 43210 

1—5 ; 

lO ' Romeel Dave 

T— I \ 

^ ' Astrophysical Sciences, Princeton University, Princeton, NJ 08544 

> ; 

^ ■ ABSTRACT 

(N ■ 
\0 ■ 

, We use hydrodynamic simulations to examine how the baryonic components of 

I galaxies are assembled, focusing on the relative importance of mergers and smooth 

I accretion in the formation of ^ systems. In our primary simulation, which mod- 

Qh' els a (50/i~^Mpc)^ comoving volume of a A-dominated cold dark matter universe, 

O ' the space density of objects at our (64-particle) baryon mass resolution threshold 

' Mc = 5.4 X W^^Mq corresponds to that of observed galaxies with L ^ -L*/4. Galaxies 

■ above this threshold gain most of their mass by accretion rather than by mergers. At 
J> . the redshift of peak mass growth, 2 ~ 2, accretion dominates over merging by about 

^ , 4:1. The mean accretion rate per galaxy declines from ~ 40Mq yr at z = 2 to 

^ I ~ IOMq yi~^ at z = 0, while the merging rate peaks later {z ~ 1) and declines more 

■ " " ' slowly, so by z = the ratio is about 2:1. We cannot distinguish truly smooth accre- 

tion from merging with objects below our mass resolution threshold, but extrapolating 
our measured mass spectrum of merging objects, dP/dM oc M~°' with a ~ 1, implies 
that sub-resolution mergers would add relatively little mass. The global star formation 
history in these simulations tracks the mass accretion rate rather than the merger rate. 
At low redshift, destruction of galaxies by mergers is approximately balanced by the 
growth of new systems, so the comoving space density of resolved galaxies stays nearly 
constant despite significant mass evolution at the galaxy-by-galaxy level. The predicted 
merger rate at z < 1 agrees with recent estimates from close pairs in the CFRS and 
CN0C2 redshift surveys. 

Subject headings: cosmology: theory — galaxies: formation 
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1. Introduction 

What galaxies look like today clearly depends on how they were assembled. For example, the 
formation of ellipticals involves considerably more random kinetic energy than does the assembly 
of spirals. Thus, it is thought that ellipticals frequently result from the rather violent merger of 
roughly equal-sized galaxies (Toomre &: Toomre 1972; Barnes & Hernquist 1992). To investigate 
this possibility, observational studies have sought to find the merger fraction as a function of 
redshift and compare this to the present abundance of elliptical galaxies (e.g., Le Fevre et al. 
1999). Conversely, the apparent quiescence of most field spirals implies tight constraints on the 
amount of merging that can occur in these systems and suggests that these galaxies more likely 
grew through a process of smooth accretion (e.g., Toth & Ostriker 1992). 

Determining galactic histories and their relation to galaxy properties is the essence of studying 
galaxy formation. Theoretical investigations of these histories must come to grips with a daunting 
array of physical processes over a wide range of scales, all of which may play important roles in 
molding a spectrum of initial density fluctuations into a population of galaxies. Hydrodynamic 
N-body simulations are a valuable tool because they can represent many of these processes and 
scales, though they still suffer from finite dynamic range and limited knowledge of the physics of star 
formation. In principle, one can use such simulations to examine the link between the formation 
histories and morphological properties of galaxies. Unfortunately, simulations that resolve the 
internal structure of galaxies and simultaneously model a representative volume of the universe 
arc (slightly) beyond the reach of current computers. Furthermore, the most ambitious efforts to 
simulate the formation of individual systems do not reproduce the observed distribution of galaxy 
properties (e.g., Navarro & Steinmetz 2000), and it is not clear whether this failure has its roots in 
incorrect cosmological assumptions, an inadequate treatment of physical processes (star formation 
and feedback in particular), or numerical limitations of the simulations themselves. 

In this paper we take a complementary approach to the study of galaxy assembly, using simula- 
tions that represent large cosmological volumes but do not resolve the internal structure of galaxies. 
Within these simulations we measure the global evolution of the number density and total mass of 
the galaxy population, and we examine the relative importance of smooth accretion and mergers 
in driving this evolution. We also investigate some related issues such as destruction of galaxies by 
tidal disruption, mass loss from galaxies during mergers, and the connection between the global star 
formation history and the history of galaxy assembly. In the context of our adopted cosmological 
model, these numerical results provide a backdrop for interpreting observations of galaxy evolution 
and merger rates or constructing theoretical descriptions of the origin of galaxy morphologies. A 
combination of N-body simulations and analytic methods has led to a fairly complete statistical 
understanding of the assembly of dark matter halos (sec, e.g., Lacey & Cole 1993; Sheth &: Tormen 
2001); our results represent an initial numerical effort to extend that understanding to the baryonic 
components of galaxies. 

This focused examination of the physics of galaxy assembly extends our earlier work on galaxy 
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formation in the framework of the inflationary cold dark matter scenario (Katz et al. 1992, 1996, 

1999; Hernquist ct al. 1995; Weinberg et al. 1997, 1999, 2000). Here we take advantage of the larger 
simulations possible with the parallel version of TreeSPH (Dave, Dubinski, &: Hernquist 1997). We 
also take advantage of an emerging consensus, driven by many independent observations, that 
favors a cosmological model with 0^ ~ 0.4, h = iifo/(100kms~^ Mpc"-*^) !=a 0.65, a fiat universe 
dominated by cold dark matter and vacuum energy, and approximately scale-invariant primeval 
fluctuations with the properties predicted by inflation. This consensus allows us to concentrate our 
efforts on a single set of cosmological parameters. However, it is worth bearing in mind that the 
viability of this model on galactic scales remains a matter of debate (see, e.g., Moore et al. 1999; 
Spergel &: Steinhardt 2000), and that even within its general framework there are uncertainties in 
parameter values that could have a significant impact on its predictions. 

The principal limitation of our study is that we cannot distinguish the process of truly smooth 
accretion from mergers with galaxies below our simulations' mass resolution threshold. Throughout 
this paper, we will use the phrase "smooth accretion" to refer to the combination of these two 
processes, as distinct from mergers with objects above our resolution threshold. For the large 
volume, high dynamic range simulation that we rely on for most of our results, the space density 
of resolved objects implies an identification with galaxies of luminosity ~ where is the 

characteristic luminosity in the Schechter (1976) luminosity function (see discussion at the end of 
§3. 

The following Section describes our numerical simulations and our method of identifying the 
dense groups of baryonic particles that represent the observable regions of galaxies. In §3 we 
investigate the global evolution of the galaxy population — number densities and mass densities 
— and we discuss numerical resolution issues in some detail. The heart of the paper is §4, where 
we examine the processes that create and destroy galaxies in the simulations, the roles of smooth 
accretion and mergers in galaxy assembly, the statistics of merger masses and mass ratios, mass 
loss during mergers, and the connection between star formation and accretion. Section 5 discusses 
the implications of these results, the comparison of our predicted merger rates to observations, and 
directions for future work. The Appendix presents details of the analysis method, which is based 
on a view of the simulations as Monte Carlo solutions of a kinetic equation describing the evolution 
of the galaxy mass function. 

2. Numerical methods 

Wc perform our simulations using the parallel version of the cosmological N-body/hydrodynamic 
code TrecSPH (Hernquist &: Katz 1989; Katz, Weinberg, &: Hernquist 1996, hereafter KWH; Dave, 
Dubinski, & Hernquist 1997), a code that unites smoothed particle hydrodynamics (SPH) (Lucy 
1977; Gingold & Monaghan 1977) with the hierarchical tree method for computing gravitational 
forces (Barnes Sz Hut 1986; Hernquist 1987). Dark matter, stars, and gas are all represented by 
particles; collisionless material is influenced only by gravity, while gas is subject to gravitational 
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forces, pressure gradients, and shocks. The gas can also cool radiatively, assuming primordial 
abundances, and through Compton cooling. In SPH, gas properties are computed by averaging 
or "smoothing" over a fixed number of neighboring particles; the calculations here use 32-particle 
smoothing. There is a maximum allowed timestep, and all particles are integrated with this step 
or one a power of two smaller. The timestep criteria are detailed further in KWH and Quinn et al. 
(1997); we set the tolerance parameter rj (defined in these papers) to 0.4. In these simulations we 
include star formation and supernova feedback using the algorithm described in KWH. For each 
star formation event, supernova energy is added to the surrounding particles as thermal energy on 
a time scale of 2 x 10^ years. 

2.1. Simulation parameters 

All four of the simulations we discuss in this paper are based on a A-dominated cold dark 
matter (CDM) cosmological model with = 0.4, Oa = 0.6, h = ifo/(100kms-^ Mpc"^) = 0.65, 
and primeval spectral index n = 0.93. With the tensor mode contribution, normalizing to COBE 
using CMBFAST (Seljak &; Zaldarriaga 1996; Zaldarriaga, Seljak, &; Bertschinger 1998) implies a 
normalization cxg = 0.8, which provides a good match to cluster abundances (White, Efstathiou, & 
Frenk 1993). We use the Hu &: Sugiyama (1996) formulation of the transfer function, and we adopt 
a baryonic density Of, = 0.02 h^"^, consistent with the deuterium abundance in high rcdshift Lyman 
limit systems (Buries & Tytler 1997, 1998), and the opacity of the high redshift Lyman-alpha forest 
(e.g. Ranch et al. 1997). All of our simulations model a triply periodic cubical volume. 

The simulations cover a range of resolutions, which allows us to examine the numerical behavior 
and test for convergence. The primary and largest simulation, L50/144 (Dave et al. 2001), has 144^ 
gas and dark matter particles in a comoving periodic box 50^~^ Mpc on a side, with a gravitational 
softening length egrav = 7 comoving kpc (equivalent Plummer softening). This simulation is 
evolved to z = 0. As discussed in §2.3 below, we adopt as our nominal gas mass resolution the 
mass corresponding to 64 SPH particles, which for this simulation is 5.4 x 10^^ Mq. We also analyze 
LI 1/64, a simulation with 64^ gas and dark matter particles in a volume ll.ll/i~^ Mpc on a side, 
also evolved to z = 0. This simulation has a particle mass eight times smaller than L50/144 and a 
gravitational softening length two times smaller. The minimum SPH smoothing length is egrav/^i 
so the higher spatial resolution extends to the hydrodynamic forces. We also analyze Lll/128, a 
simulation with 128'^ particles of each type in a volume the same size as the Lll/64 simulation, but 
only evolved to z = 3. It has a particle mass eight times smaller than Lll/64 and a gravitational 
softening length two times smaller. Finally, we analyze Lll/64'. This simulation has the same 
resolution and volume as Lll/64, but it uses the same random realization of the CDM power 
spectrum as Lll/128. It has been evolved to z = 3, and the comparison to Lll/128 gives our most 
direct test of finite resolution effects down to this redshift. 

Our results depend on the spatial and mass resolution of the simulation, which we discuss in 
more detail below, and on the temporal spacing of the available outputs. Large time intervals tend 
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to increase the amount of smooth accretion because small groups can form and merge with larger 
neighbors within a single time interval. For our analysis, we chose a sequence of outputs with a 
spacing of Az = 0.5 for 7 > z > 4, Az = 0.25 for 4 > z > 1, and Az = 0.125 for 1 > z > 0. 
The corresponding time intervals, as seen in the figures below, are roughly 1/3 Gyr, 2/3 Gyr and 
4/3 Gyr, respectively. These are on the order of the infall time scale at the corresponding epoch, 
although perhaps somewhat longer at the latest time. We repeated some of our main analyses using 
twice as many outputs, i.e. halving the redshift spacing, and found similar results. 

2.2. Group finding algorithm 

SPH simulations with CDM initial conditions and radiative cooling lead to the formation of 
dense groups of baryonic particles that have sizes and masses comparable to the luminous regions of 
observed galaxies (Katz et al. 1992; Evrard, Summers, & Davis 1994). If star formation is included, 
these dense groups are the regions where stars form (KWH). The accretion and merger histories of 
these objects are the subject of this paper, so the identification of distinct particle groups underlies 
all of our subsequent analysis. We use the group finding algorithm of Stadel et al. (2001), Spline 
Kernel Interpolative DENMAX, abbreviated SKID (Gelb & Bertschinger 1994; see also KWH and 
http:/ /www-hpcc. astro. Washington. cdu/TSEGA/tools/skid. html). The basic algorithm consists of: 
1) determining the smoothed density field; 2) moving particles upward along the gradient of the 
density field using an heuristic equation of motion that forces them to collect at local density 
maxima; 3) defining the approximate group to be the set of particles that aggregate at a particular 
density peak; 4) finally, removing particles from the group that do not satisfy a negative energy 
binding criterion relative to the group's center of mass.^ In contrast to the identification of dark 
matter halos in N-body simulations, there is essentially no ambiguity about the identification of 



'^This procedure starts with the highest energy particles and updates the potential as particles are lost. 
Table 1: Simulations 



Name 


L(Mpc) 


N 




MmaxiMQT 






L50/144 


50/i-i 


2 X 144^ 





3.1 X 10^2 


6.8 X 10^ 


5.4 X 10^° 


Lll/64 


llh-^ 


2 X 64^ 





3.2 X 10^° 


8.5 X 10^ 


6.8 X lO'^ 


Lll/128 


llh-^ 


2 X 128^ 


3 


3.0 X 10^° 


1.1 X 10^ 


8.5 X 10^ 


Lll/64' 




2 X 64^ 


3 


3.0 X 10^° 


8.5 X 10^ 


6.8 X lO'^ 



"Mass of lO"' largest galaxy in simulation, at z 
''Mass of 8-particle group 
■^Mass of 64-particle group 



= for L50/144 and a = 3 for other simulations 
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distinct baryon clumps in these SPH simulations because dissipation greatly increases the density 
of cooled baryons with respect to the local background. 

As output, SKID produces a list of baryonic particles, both gas and stars, that belong to each 
identified group. To be included in a group, gas particles must have a temperature T < 3 x lO^K 
and and a density p > lO^Jlfe. In effect, we consider only the cold gas and stars, the material that 
comprises the bulk of the matter in the central, visible regions of galaxies. Therefore, throughout 
this paper we use the terms group and galaxy interchangeably. 

One complication that arises with this group definition is lost particles. Specifically, given 
all the particles in groups at one output, SKID finds that roughly 1% of these are lost by the 
next output. However, more careful inspection shows that roughly 30% of those particles declared 
lost still reside physically within the group boundary, which we define to be the maximum radius 
of any bound member particle. These particles may in fact be unbound and in the process of 
leaving the group, although we have not checked their energies relative to the group center-of- 
mass. Nevertheless, in the interest of estimating the mass loss conservatively, we simply add these 
particles back to the group. Consequently, our group definitions become dependent on output 
spacing, although the effect is only of order 1%. 

2.3. Mass resolution 

We allow SKID to identify groups with as few as eight particles. Because these particles must 
be stars or cooled baryons, even these small groups are real, high contrast objects residing in dense 
backgrounds. However, our comparisons, statistical and direct, between simulations with different 
mass resolutions show that a simulation's group list becomes incomplete (relative to a simulation 
with higher resolution) below ~ 60 particles. Above this threshold, the total baryon mass of groups 
is well converged (see the discussion of Figure 1 below), though the division into stellar and gas 
components remains resolution dependent up to a somewhat higher threshold. For most of the 
analyses of this paper, we adopt Ngrp = 64 particles as the minimum number for a resolved group, 
and we discard all SKID groups with fewer particles. Our nominal mass resolution limit is therefore 
that of 64 SPH particles, though since mass is redistributed among particles during star formation 
(see KWH) the exact mass of a 64-particle group may be slightly higher or lower. We further 
discuss resolution effects and other limits on the mass ranges probed by our simulations in the next 
Section. 

3. Evolution of galaxy number and mass densities 

First we look at the simplest quantities that characterize the galaxy populations in the sim- 
ulations: the mass function, the total number density, and the total mass density. Examining 
these characteristics shows qualitatively how the galaxy population evolves and also illustrates the 
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numerical factors that come into play in the analysis. 

As mentioned above, the most important numerical variable is the mass resolution, the min- 
imum mass of a resolved particle group. Because the number of galaxies increases towards low 
masses, the number of galaxies in a simulation depends sensitively on this lower mass limit. The 
box size also comes into play, however, as this determines the number of relatively rare, high- 
mass objects that form, thereby determining the maximum mass for which statistical results are 
meaningful. 

Table 1 lists the numerical parameters of the simulations described in §2.1, in particular the 
characteristic minimum and maximum masses associated with the resolution and box size. We 
define the maximum mass, Mjnax, to be the baryonic mass of the 10*'* largest galaxy in the simula- 
tion volume, identified as described in §2.2. M^nin denotes the mass corresponding to an 8-particle 
minimum group size, and Mc denotes the adopted cutoff mass corresponding to the 64-particle 
minimum. All three ll.ll/i~^ volumes have about the same maximum mass at z = 3. The min- 
imum mass is determined by the resolution, so the Lll/128 simulation has a substantially lower 
minimum mass. 

Figure 1 shows particle groups in the L50/144, LI 1/64, and Lll/128 simulations at z = 3. 
In any given simulation, the mass function begins to turn over at a mass 5-10 times higher than 
the 8-particle minimum of the SKID groups. Comparison to the next higher resolution simulation 
shows that this turnover is a numerical artifact. If we instead restrict attention to masses above our 
adopted 64-particle threshold Mc, then the mass function in a higher resolution simulation gives 
a fairly smooth continuation of that in the lower resolution simulation. The envelope of the three 
mass functions approximately follows a Schcchter function. 

All the analyses that follow include only galaxies with masses above Mc- Because of their 
different resolutions and box sizes, our simulations probe different mass ranges. In the case of 
L50/144, one can see from Figure 1 that this range is nearly disjoint from that of the smaller 
volume simulations. 

The lower panel of Figure 2 compares the comoving number density f of galaxies above the 64- 
particle threshold in the four simulations, as a function of time (lower axis label) or redshift (upper 
axis label). The upper panel shows the comoving density /j, of baryonic mass in these resolved 
galaxies. Once again, the effects of mass resolution and box size are clearly evident. Lll/128 has 
the highest number and mass densities. LI 1/64 and LI 1/64' have fewer galaxies, of course, since 
they have the same upper mass cutoff as Lll/128 but a lower mass cutoff that is eight times higher. 
However, if the power law form of the low end of the mass function, roughly n(M) cx M~^'^ in 
Figure 1, continued to arbitrarily high mass, then L50/144 would have the highest number and 
mass densities because it has the largest dynamic range Mmax/Mc- (For n oc M~^, the number 
of galaxies per logarithmic interval is constant, and the mass per logarithmic interval increases 
linearly with M.) In reality, the maximum mass M^ax of L50/144 lies well into the exponential 
cutoff regime of the mass function, and as a consequence the galaxy population of L50/144 has the 




Fig. 1. — The mass function at 2: = 3 in the L50/144 (sohd histogram), LI 1/64 (dashed histogram), 
and Lll/128 simulations (dotted histogram). Error bars reflect Poisson errors. The solid curve is 
a Schechter function overlaid for comparison. 
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lowest mass density and by far the lowest number density of the four simulations. 

Turning from numerical considerations to evolution, we see that number and mass densities 
both increase most strongly at early times, i.e. for z > 2. In the simulations that continue to 
z = 0, the number of objects becomes approximately constant for z < 2, and the difference in 
mass density between L50/144 and Lll/64 stays about the same. As discussed in §4 below, the 
number of galaxies remains roughly constant even though merging continues to late times because 
new galaxies also continue to form. In contrast to the number density, the total mass in galaxies 
continues to grow, both from the formation of new groups and from the accretion fed growth of 
existing galaxies. 

Figure 3 shows this behavior in terms of the corresponding rates of change in mass density (i 
and number density z>. The differences between simulations are greatest at high redshift, where 
the hierarchical mass scale is lowest. For z < 2 the results at different resolutions agree, indicating 
that the formation rate of very low mass objects has become negligible. The differences in total 
mass and number density apparent in Figure 2 must therefore be imprinted at early times. 

The physical quantity one would like to know is the mass or number density of galaxies above 
some specified mass M. In practice AI cannot be below the mass resolution of the simulation, 
Mc, and the mass and number densities measured by the simulations are not those for all galaxies 
above a given mass but only for galaxies below a maximum mass determined by the volume. 
However, we can measure the same quantity in more than one of our simulations to assess the 
robustness of our measurements given these limitations. The Lll/64 and Lll/64' simulations 
have the same resolution and box size, but there is a small difference in the quantities plotted in 
Figures 2 and 3, represented by the filled squares and four point stars. These differences show 
the random fluctuations associated with two different realizations of the initial power spectrum in 
this volume. We can assess the impact of our finite resolution on these quantities by comparing 
the Lll/64' simulation (four point stars) to the LI 1/128 simulation analyzed with the mass cutoff 
Mc = 6.8 X 10^ Mq of Lll/64'. These results are shown by the three point stars in Figures 2 and 3. 
The good agreement between the three and four point stars indicates that we can robustly measure 
number and mass densities for galaxies above the 64-particle threshold mass Mc, and that finite 
resolution effects on our statistics above this threshold should be small. 

To relate our numerical results to an observational context, we would like to know the ap- 
proximate luminosity of galaxies corresponding to our resolution threshold Mc- Although we have 
stellar masses for each of our galaxies, the mass-to-light ratio of the stellar populations depends 
sensitively on the assumed initial mass function and on the age, metallicity, and extinction of the 
stellar population. Furthermore, the overall scale of the galaxy baryon masses may be sensitive 
to some cosmological parameters, especially O5. Given these uncertainties, the most robust way 
to identify our mass threshold with an approximate luminosity threshold is by matching the space 
densities of simulated and observed galaxy populations. If we match the mass function in Figure 1 
to Blanton et al.'s (2001) determination of the luminosity function from the Sloan Digital Sky Sur- 
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Fig. 2. — The total comoving mass (top) and number density (bottom) of galaxies as a function 
of cosmic time for galaxies with masses above Mc = 5.4 x 10^° Mq (L50/144-open squares), = 
6.8 X IO^Mq (Lll/64-solid squares, Lll/64'-four point stars, Lll/128-three point stars), and 
Mc = 8.5 X 10^ Mo (Lll/128-open triangles). 
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Fig. 3. — The rates of change of mass (top) and number density (bottom) above the masses shown 
in Figure 2. The number and mass densities for masses above = 5.4 x IO^^Mq (L50/144-open 
squares) and = 6.8 x 10^ (Lll/64-soUd squares) agree for z < 1.5. 
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vey, the implied luminosity corresponding to the threshold Mc of L50/144 is ~ -L*/4, the precise 
value depending on whether we match the g' or r' luminosity function and whether we match space 
densities at Mc or at L*. 

4. Evolution of the galaxy population through accretion and mergers 

4.1. Galaxy creation and destruction 

The instantaneous rate of change in galaxy number in the simulations is determined by the 
relative rates of formation and destruction. Galaxies are destroyed either through merging with 
another larger galaxy or though disruption. In the upper panel of Figure 4, open squares show 
the rate of change z> of galaxy number density from the L50/144 simulation. (These are the same 
as the open squares in the lower panel of Figure 3, but the vertical scale is greatly expanded.) 
Filled squares show the creation rate of new galaxies; in general these galaxies are not entirely 
new, but they have gained enough mass since the last output to cross the mass threshold Mc- The 
difference between the filled and open points is the destruction rate, shown also by the filled squares 
in the lower panel. Open squares and filled triangles in this panel show the contributions to the 
destruction rate from mergers (which make one galaxy from two) and disruption (which moves a 
galaxy from above Mc to below Mc), respectively. Figure 5 shows creation and destruction rates 
for Mc = 6.8 X 10^ Mq, from Lll/64. This plot is noisier because of the smaller number of galaxies 
in the simulation, but it is qualitatively similar except for the higher number densities associated 
with less massive galaxies (note the change in vertical scale from Figure 4). 

In both cases the creation rate climbs rapidly to a peak at z ^ 3 and declines steadily thereafter, 
though new galaxies continue to form (i.e., to cross the Mc threshold) down to z = 0. The merger 
rate climbs more gently, reaching a broad maximum at z 1 — 2 and declining only slowly at lower 
redshift. Mergers always dominate over disruption as a destruction mechanism. At high redshift 
the creation rate is much larger than the destruction rate, but the balance begins to shift at z < 3. 
By the present day, the creation and destruction rates are nearly equal, and both are much larger 
than the net rate. This result shows that a high galaxy merger rate need not lead to rapid evolution 
of the galaxy luminosity function, since other galaxies can grow to replace the ones that are lost. 

4.2. Growth of mass through accretion and mergers 

Figures 2 and 3 show that the total mass in resolved galaxies increases with time, as one would 
expect. Since we only count mergers among resolved objects, all new mass must enter the population 
through accretion. However, mergers redistribute mass within the resolved galaxy population, and 
we can sensibly ask whether existing resolved galaxies gain more of their mass through accretion or 
through merging with other resolved systems. Our methods for determining the smooth accretion 
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Fig. 4. — Contributions to the number density evolution for galaxies with masses above Mc = 
5.4 X 1O^°M0 (L50/144). Top panel: filled squares give the creation rate of new galaxies, and open 
squares give the net rate (creation minus destruction). Bottom panel: filled squares give the total 
destruction rate, and open squares and filled triangles show the respective contributions of mergers 
and disruption. 
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and merger rates from the simulation outputs are fully described in the Appendix. The most 
important point to recall here is that we define accretion to be true smooth accretion plus any 
merging with galaxies below the cutoff mass Mc, since we have no reliable way to distinguish these 
processes. For the remainder of this paper we will restrict our discussion to the L50/144 simulation 
because it has the largest dynamic range and probes the mass range of greatest interest. 

In Figure 6, pentagons show the total mass accretion rate for galaxies with masses above 
Mc = 5.4 X IO^^Mq; except for the minor effects of disruption and mass loss, this is the same 
quantity shown by the open squares in Figure 3 (top panel). Since we would like to know how 
existing galaxies acquire mass, we subtract the contribution from newly formed galaxies and plot 
the remainder as the filled squares. This globally averaged mass accretion rate in galaxies above 
5.4 X 10^^ Mq rises rapidly until z ~ 2 and declines slowly thereafter. At z = 0, the accretion rate is 
about 1/3 of its peak value. The rate of truly smooth accretion by galaxies above Mc is necessarily 
lower than that shown in Figure 6, since higher resolution simulations would resolve some of this 
accretion into mergers with small groups. In addition to the total accretion rate, we show the 
contributions from accreted gas (open squares) and stars (crosses) separately. As expected, gas 
dominates the type of accreted material. 

Figure 7 shows the distribution of accretion rates at four different redshifts. The mean mass 
accretion rate of resolved galaxies drops from M ~ 4OM0yr~^ atz = l — 2toM~ lOM0yr~"^ at 
z = 0, though at each redshift the M distribution is broad. Gas accretion always dominates over 
stellar accretion, as expected from the global properties in Figure 6. 

Figure 8 shows the volume averaged rate at which resolved galaxies (M > 5.4 x IO^^Mq in 
L50/144) gain mass through merging with other resolved galaxies. The single most important 
result of this paper comes from the comparison of Figure 8 to Figure 6: galaxies gain most of their 
mass by accretion, not by mergers (note the large change in vertical scale). At 2; ^ 2, accretion 
dominates merging by about a factor of five. However, accretion declines more rapidly than merging 
towards low redshift, and by 2; = it dominates by only a factor of two. The ratio of accretion 
growth to merger growth in galaxies above this mass threshold would drop if our simulations had 
higher mass resolution and could therefore resolve objects that are currently counted in the smooth 
accretion rate.^ Nevertheless, the total merging rate is fully determined for galaxies with a mass 
above Mc = 5.4 x lO^^M©, and in this mass range mergers contribute < 1/3 of the mass growth 
rate at every redshift. Furthermore, we will show in §4.3 that most of the merger contribution 
comes from relatively massive objects, so the overall accretion-to-merger ratio would remain high 
even with higher mass resolution (see discussion in §5). 

Focusing on the merger rate itself, we see that it climbs fairly quickly until z ~ 1.5. The 
variations thereafter appear largely stochastic, with a slight overall decline. The clear decline in 



^However, a higher resolution simulation would not necessarily give a lower accretion-to-merger ratio for its full 
population of resolved galaxies, since its mass resolution threshold would also be lower, and these less massive galaxies 
would still tend to accrete unresolved material. 
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Fig. 6. — Evolution of the mass accretion rate for galaxies above Mc = 5.4 x lO^^M© (L50/144). 
solid pentagons show the total including the contribution from newly formed galaxies. Solid squares 
show only the accretion rate onto existing galaxies, excluding newly formed objects. Open squares 
show the contribution to this accretion from gas and crosses the contribution from stars. 
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Fig. 7. — Distribution of accretion rates onto individual galaxies at different rcdshifts for galaxies 
■with masses above Mg = 5.4 x IQ^^Mq (L50/144). The solid line shows total accretion rates, the 
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Fig. 8. — The rate of mass gain through merging as a function of time, for galaxies with masses above 
Mc = 5.4 X 10^° M0 (L50/144). Sohd squares show the total rate, open squares the contribution 
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factor of five compared to Figure 7. 
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the number rate of merging (Figure 4) is compensated by an increase in the typical mass of the 
merging objects. Open squares and crosses in Figure 8 show the separate contributions of gas and 
stars. Gas-rich mergers are important early on, where about half the accreted material is gas and 
half is stars, but they are quickly overtaken by mergers with predominantly stellar systems. 

4.3. Masses and mass ratios of merging galsLxies 

Figure 9 shows the probability distribution of the masses of smaller "satellite" galaxies that 
merge into larger "parent" galaxies. Above the cutoff mass, M^, the distribution generally appears 
to follow a power law in mass. The slope of this power law decreases with redshift, presumably 
following the evolution of the non-linear mass scale. At high redshift {z > 0.5), the estimated slope 
is steep enough (a < — 1) that small objects dominate in number, while large objects dominate in 
mass (i.e., an integral for the total mass of merging objects diverges at the upper end). However, 
at low redshift {z < 0.5), the slope flattens, so that large objects dominate in both number and 
mass. 

Figure 10 shows the distribution of merger mass ratio, the ratio / = Mgat/^par between the 
satellite galaxy and the larger parent with which it merges. At the earliest times {z > 2), large mass 
ratios predominate, although the total amount of merging is somewhat smaller than at later times. 
Subsequently, a low / tail develops, and the distribution becomes approximately scale free above the 
effective resolution cutoff. The probabilities extend to lower and lower / as the simulation evolves. 
The slopes of the distributions arc quite steep, with higher mass ratios (/ > 0.4) dominating by 
number. However, it is important to keep in mind that our finite resolution suppresses low-/ 
mergers, and since the number of resolved galaxies is largest near Mc, there is some preference for 
merger ratio near unity almost by definition. 

Given the resolution-limited mass, Mc, the simulation also gives an exact or complete estimate 
of the total amount of merging above a given mass ratio and corresponding parent mass. Stated 
in more mathematical terms, for a given merger ratio /o, there is a corresponding parent mass, 

Mq = Mc/ fo- for which the simulation includes all mergers with / > /o and Mpar > Mq. Conversely, 
given /o, the rate is not exact for Mpar < Mq because the total should include satellites with 
Msat < Mc, which are unresolved in the simulation. 

Figure 11 shows the total number merger rate t'(> Mq, > Jq) at different redshifts, while Figure 
12 shows the corresponding mass merger rate /i(> Mq, > /o). Both figures display a range of values 
of /o and Mq. Of particular interest is the amount of major merging, since that could transform 
a galaxy's morphological type. We define a major merger to be / > 0.25.^ With the present 



^Simulations show that for / ~ 1, mergers essentiaUy transform galaxy types (e.g. Barnes 1988, 1992; Hernquist 
1992, 1993; Barnes & Hernquist 1996), while for / ~ 0.1, galaxies are damaged, but not completely destroyed (e.g. 
FuUager et al. 1993; Walker et al. 1996). 
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Fig. 9. — The mass distribution of galaxies that merge into larger galaxies as a function of smaller 
galaxy mass, in four redshift ranges, for galaxies with masses above Mc = 5.4 x IO^^Mq (L50/144). 
The solid line shows total merged mass, the dotted line the gas mass, and the dashed line the stellar 
mass. There is a cutoff in total mass at logM = 10.8, corresponding to the minimum mass Mc, 
although stellar and gas contributions can be smaller. The line indicates the approximate slope a 
of the high-mass end. The arrow on the abscissa denotes the mean merger mass. 
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simulation, we can measure the amount of major merging exactly for Mpar > Mq = 2.2 x 10 Mq, 
approximately L*. Overall, major merging appears to contribute significantly to the total amount 
of merging measured in the simulation. For this threshold mass at z ~ 1, 50% of the total mass 
in merging comes from pairs with /o > 0.25. At z = 0, the contribution increases to nearly 75%. 
In terms of number, the rate lies in the range from 10~^ — 10~^ Gyr~^ Mpc~^ for z < 3. Figure 13 
shows the rates for both number and mass as a function of redshift for major merging. For z > 3, 
major merging does not occur at this resolution limit. For z < 3, both rates climb. The number 
rate reaches a stochastically varying plateau for z < 2, while the mass rate continues to climb until 
z = 1, where it also reaches a plateau. 

4.4. Mass loss 

Fairly significant mass loss also occurs in the simulations. Figure 14 shows the mass lost from 
galaxies as a function of time for galaxies with masses above = 5.4 x IO^'^Mq (L50/144). The 
trend follows the merging history shown above — it rises until z ~ 1 and reaches a plateau for z < 1. 
Overall, mass loss rates are about 25% of merging rates. Gaseous mass loss is nearly constant with 
time, while stellar mass loss increases with time. Stellar mass loss dominates after z ~ 0.5. Visual 
inspection shows that the mass loss occurs in highly clustered regions. 

4.5. Star formation 

The volume- averaged star-formation rate is closely related to the growth of galaxies in the 
simulation. Figure 15 plots the global star-formation rate figf for galaxies with masses above 
Mc = 5.4 X 10^^ Mq (L50/144) as a function of time. Interestingly, the overall rate (filled squares) 
closely follows the global rate of mass accretion shown in Figure 6 rather than the merging history 
shown in Figure 8. We also plot separately the star formation contributed by gas that is newly 
acquired through smooth accretion since the last time output and that contributed by gas already 
present within the parent galaxy or merged satellite in our last time slice (crosses and open squares, 
respectively). Most of the stars form from gas already present in galaxies. Moreover, since merging 
rates are much smaller than star-formation rates (c.f. Figures 6 and 8), this implies that star 
formation typically occurs in gas that was present in the parent galaxy at the previous timestep. 
Hence, in our simulations, gas that enters a galaxy typically waits for at least the interval between 
our outputs (1/3 Gyr at high redshift and 4/3 Gyr at low redshift) before converting to stars. 

5. Discussion 

The results presented here describe the growth of the central baryonic components of galaxies 
in the context of hierarchical structure formation. We have focused primarily on the relative 
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Fig. 11. — The cumulative number merger rate as a function of parent group mass for different 
merger ratios: i'{> Mq, > /o). Each sequence of points indicates a different fo- from top to 
bottom, fo = 1%, 25%, 50%, 75%, 100%. The solid symbols indicate the parent mass range for 
which Msat > Mc for / > fo, i.e. the mass of the merged object lies above the 64-particle threshold. 
The open symbols indicate parent masses for which Mgat < Mc- Alternating symbols are used for 
visual clarity. 
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Fig. 14. — The rate of mass loss from existing groups per unit comoving volume as a function of 
time, for galaxies above = 5.4 x lO^^M© (L50/144). Solid squares show the total, open squares 
the contribution from gas, and crosses the contribution from stars. 
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Fig. 15. — The star formation rate in galaxies above = 5.4 x 10^^ Mq (L50/144) per unit volume 
as a function of time. The solid squares show the total, the open squares the contribution from gas 
already in groups and the crosses the contribution from gas that is recently accreted. 
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contributions of merging and smooth accretion to the rate at which large galaxies gain mass. 
Overall we find that accretion dominates, especially at higher redshift, while merging becomes 
progressively more significant at later times. 

The principal limitation of our results is that wc cannot distinguish between smooth accretion 
and merging with objects below our mass resolution threshold, which is Mc = 5.4 x 10"'^'' Mq for 
our primary simulation (L50/144). As discussed at the end of §3, matching the observed galaxy 
space density implies that this mass threshold corresponds approximately to a luminosity L*/4. 
Our analysis gives well defined, and, we believe, numerically robust predictions for the merger rates 
of galaxies above this threshold with other galaxies above this threshold. These are lower limits 
to the merger rates of galaxies above this threshold with all other galaxies, and our corresponding 
accretion rates are upper limits to the rates of truly smooth accretion by these galaxies. 

The mass spectrum of merging galaxies shown in Figure 9 gives useful guidance to the pos- 
sible contribution of sub-resolution merging. Over roughly a decade in mass, the mass spectrum 
is approximately a power law oc M", with a ~ —0.8 at z = and a ~ —1.4 at z = 1 — 3. The 
cumulative distribution of mergers above mass M scales as M"+-'^, and the amount of mass con- 
tributed by mergers above mass M scales as M"+^, thus diverging towards the high mass end of 
the spectrum. Thus, the measured slopes of the merger mass spectrum indicate that sub-resolution 
mergers should contribute relatively little mass compared to resolved mergers. To obtain a quanti- 
tative estimate, we assume that the measured power law mass spectrum extends from a high mass 
cutofi" Mh = 6 X W^^Mq (see Figure 9) down to a low mass cutoff M; = 0. We then find that 
merging galaxies below Mc should contribute roughly 25% of the total mass in mergers for 2; > 1 
(a ^ —1.4) and only about 5% of the total for z < 0.5 {a ~ —0.8). Thus, at high redshift, the 
simulation appears to account for more than 75% of the total mass in merging, whereas at low 
redshift it accounts for at least 95%. Combined with the results in Figures 6 and 8, this implies 
that merging overall accounts for no more than 25% of the total mass in accretion and merging 
combined for z > 1 and no more than 35% of the total for z < 0.5. While higher resolution sim- 
ulations will be needed to verify this estimate, the extrapolation from our current results implies 
that truly smooth accretion always exceeds merging by at least a factor of two in contributing to 
the mass evolution of galaxies. 

It is interesting to compare the mass spectrum of mergers of baryonic objects found here to the 
mass spectrum of dark matter halo mergers, derived by Lacey & Cole (1993) using the extended 

Press-Schechter formalism. The comparison is not exact because we determine the mass distribution 
of all merging objects, whereas Lacey &: Cole (1993) determine the distribution of objects falling 
into a single halo. The merger mass spectrum that they find is steeper than the one shown in our 
Figure 9. Nevertheless, the high redshift behavior is qualitatively similar in both cases: low mass 
objects dominate in number but high mass objects dominate in mass. At low redshift there is a 
qualitative difference in the two rcsTilts: our scaling implies that high mass objects dominate in both 
mass and number, while Lacey & Cole's implies that low mass objects still dominate in number but 
high mass objects dominate in mass. This difference could plausibly refiect the different dynamics 
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of baryon and dark matter mergers. When two dark halos merge, the galaxies that they contain 
will not merge immediately. Dynamical friction can drag together massive galaxies more quickly 
than low mass galaxies, and this should make the merging mass spectrum of galaxies shallower 
than that of dark halos. 

Our main point of contact with observations is through merging rates. Observationally, the 
efforts focus on determining the evolution of pair fraction with redshift (Zepf &: Koo 1989; Carlberg 
et al. 1994; Yee & Ellingson 1995; Woods, Fahlman, & Richer 1995; Patton et al. 1997; Abraham 
1998; LeFevre et al. 1999; Carlberg et al. 2000). The difficulty lies in identifying pairs of galaxies 
that should merge in a time of order the system dynamical time at increasingly large redshift. Here 
we use the true physical merger rates, since these are the quantities we have available from the 
simulations, but we do not investigate whether the observational analyses accurately infer these 
rates. 

From their analysis of pairs selected from the CN0C2 and CFGRS surveys, Carlberg et al. 
(2000) have recently estimated the rate of mass growth due to merging in galaxies above about 
0.2L^, out to z = 1. They find an integrated mass growth rate of about 2 x 10~^L* Gyr~^ per 
galaxy, with an estimated uncertainty of a factor of two. Their 0.2L* threshold should correspond 
approximately to our Mc threshold, though there is a factor ~ 2 uncertainty in this identification. 
While Carlberg et al. (2000) do not explicitly present mass ratios for their pairs, it is likely that for 
such large galaxies the mergers are predominantly major (/ > 0.25). From Figure 13 we find that 
the major merging rate fimge at and above varies from 2 x 10~^ to 3.5 x 1O~^M0 Mpc~^ yr~^ 
between z = 1 and z = 0, while the comoving number density of objects varies from 4.2 x 10~^ to 
5.5 X 10~^ Mpc~^ over the same interval. Dividing the merger rate by the comoving number density 
and by an assumed value of 2.2 x 10^^ MqL^^ yields rates of mass growth from major merging in 
the range 2 x 10~^ — 5 x lO^^L* Gyr^^ per galaxy. The agreement with the Carlberg et al. (2000) 
results is encouragingly good, given the theoretical and observational uncertainties that enter the 
comparison. 

Another interesting feature of our results is the mass loss due to merging. Although the group 
finding algorithm introduces some uncertainty in the measurement, it does appear that a significant 
fraction of the mass in merging satellites winds up as intergalactic debris. Visual inspection of the 
simulations shows that this material permeates galaxy clusters and large galaxy groups and forms 
low-density halos around massive galaxies. Recent observational and theoretical work has provided 
evidence for the existence of such material, particularly in galaxy clusters (e.g. Thcuns &: Warren 
1997; Ciardullo et al. 1998; Ferguson, Tanvir & von Hippcl 1998; Calcanco-Roldan et al. 2000). 
In future work, we will examine in more detail the properties of the debris produced by merging 
galaxies. 

One obvious way to extend the work presented here is to employ new, higher resolution sim- 
ulations. We are currently evolving a simulation similar to Lll/64 but using 128^ particles in 
a 22.22^"^ Mpc box. This simulation will have larger dynamic range than Lll/64 and greater 
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overlap with L50/144, improving our ability to assess numerical resolution effects. In the longer 

term, a simulation with 288^ particles in a 50h~^ Mpc box would have the same mass resolution, 
Mc = 6.8 X IO^Mq, as Lll/64, sufficient to settle many of the remaining issues regarding the growth 
of galaxies in numerical simulations. 

Motivated in part by the kinetic equation approach described in the Appendix, we have focused 
here on volume-averaged rates of accretion and mergers. These rates provide a good global view of 
activity in the resolved galaxy population, but there are other questions that can only be addressed 
by examining accretion and merger histories as a function of environment or galaxy by galaxy. 
For example, we would like to know what fraction of galaxies experience a major merger between 
z = 0.5 and the present day, or what fraction have quiescent accretion histories from z = 1 to 
z = 0, and we would like to know how those fractions correlate with galaxy mass, environment, and 
stellar age. The analysis of mass acquisition can also be extended to examine how galaxies gain 
angular momentum or random kinetic energy. We will present results from such analyses in future 
work. Studies of this sort, applied to steadily improving numerical simulations, should provide a 
solid foundation for understanding the mass, luminosity, and morphological evolution of the galaxy 
population. 

Wc thank Jeff Gardner for providing the basis for our merger code and for performing the 
LI 1/128 simulation. We also thank Mark Fardal, Enrico Vesperini and Martin Weinberg for helpful 
discussions. This work was supported by NASA Astrophysical Theory Grants NAG5-3922, NAG5- 
3820, and NAG5-3111, by NASA Long-Term Space Astrophysics Grant NAG5-3525, and by the 
NSF under grants ASC93-18185, ACI96-19019, and AST-9802568. The simulations were performed 
at the San Diego Supercomputer Center and NCSA. 



A. Evolutionary equation 

The simulations analyzed above can be viewed as Monte Carlo sohitions of the kinetic equation 
that describes the evolution of the mass function n{M, t) (and includes a large number of additional 
degrees of freedom, which we have projected over in our analysis). In the continuum limit, a simple 
form of the equation can be written: 

— = - — \nM]+ / dM'n{M - M',t)n(M',t)r(M - M',M',t) 
ot dM Jm^ 

/>oo 

- n{M,t) / dM'n{M',t)T{M,M',t) + C{M,t)-D{M,t). (Al) 

On the right-hand side, the first term gives the rate of change of the baryonic mass M of a galaxy 
with accretion rate M. Accretion gives rise to advective evolution in the mass phase space so the 
term is analogous to a continuity term. In writing this, we have assumed that the accretion M 
depends uniquely on the mass M, whereas, in reality, a sample of galaxies at mass M would have 
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a distribution of accretion rates. What is written here can be thought of as the evolution of the 
'average' galaxy of mass M. 

The second and third terms describe the merging of galaxies. The former gives the rate of 
creation of new galaxies of mass M through the merging of pairs with masses M' and M — M' at 
the rate T{M — M',M,t). The latter gives the loss rate of galaxies at M through collisions with 
galaxies at M' at the rate r(M,M',t). Note that we have ignored mass loss and disruption in 
writing these terms. 

The fourth and fifth terms denote the creation and destruction of galaxies of mass M, respec- 
tively. In the simulation, galaxies are not created at all masses; instead they usually pass directly 
through the minimum mass threshold Mc by mass accretion. In this case, one could introduce 
the creation rate as a boundary condition. Sometimes, however, sub-Mc galaxies can merge and 
produce a new galaxy with a mass somewhat larger than Mc- As a result, we can qualitatively 
describe the creation function as a onc-tail distribution that peaks at M,. and has a fairly short 
tail to higher mass. Pure destruction does not really occur in the simulations at our adopted mass 
threshold, and the term has only been included for completeness. 

By taking the first moment of equation (Al) with respect to M, one obtains the equation 
describing the evolution of the total baryonic mass in galaxies. From this, one can see that the 
total mass can change only through accretion onto and mass loss from existing galaxies and through 
creation and destruction. Merging, by contrast, alters only the number of galaxies. 



It is straightforward to generalize equation (Al) to the discrete form that is required for analysis 
of the simulations. Assuming the discretized form of the equations, we write down the expressions 
relating the total accretion and mass loss to the change in the total amount of mass in galaxies 
between two times. Let there be A^i groups at time ti and N2 groups at time t2 where t2 > h. 
Then the total mass in groups at either time tk 



B. Changes in total mass 




(Bl) 



I 



and the change in mass in groups between the two times 



AM = M2 - Ml = Mace - Ml, 



(B2) 



where Mace is the mass accreted smoothly and Mi^ss is the mass lost from the progenitor groups 
at time ti. Here we have included mass growth from created galaxies in the accretion term and 
mass loss from destroyed galaxies in the loss term. Defining /j as the fraction of group i at ti that 
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contributes to some group at t2, we can write 



Mioss = 51(1 - fi)mi, (B3) 

i 

So that 

Ni 

Mace = M2 - (Ml - Mioss) =M2-J2 fi^i = M2 - Mgrp, (B4) 

i 

where Mgrp denotes the group mass at ti that contributes to the group mass at t2- Using Mg^p, 
we can define the merged mass 

N2 

Mmerge = Mgrp — ^ ^ fimaxj''^imaxj-> (■^^) 
j 

where fi^^^j and rrii^^^j denote the fraction and total mass of the largest group at t\ that contributes 
to group j at t2- 



C. Changes in gas/star mass 

The expressions given above change when considering individual gas and star components 
because of star formation. In this case, defining Mj as the total gas mass in iVj groups at time 
we can write the total change in gas mass in groups: 

AM^ = M| - Mf = Mf,, + Ml - Ml^^ - (CI) 

where Mace denotes the amount of gas accreted into groups, Mf the gas mass created from stars, 
^loss amount of gas lost from pre-existing groups, and Mg the stellar mass created from gas. 
Of course, in the present simulations there is no stellar evolution mass loss so M| = 0. As in the 
definition above, Mf^^^ = 5^j(l — fi)'mf, the respective sum of the fraction of mass in groups at 
ti that does not remain in groups at t2 ■ For stars, the treatment is analogous; however the star 
formation terms change sign. 

To determine the Mg term, we must consider the gas and star masses of the individual particles 

at consecutive outputs. However, because of the finite time resolution, it is impossible to say 
whether star formation occurs before or after material is added to the new system. Here we assume 
that star formation occurs after material is added to a new system at time ^2- 

With the total particle mass at time ti fik,i = j-i-f. i{ti) + /x^ i, the total stellar mass in groups 
at ti that contribute to the total mass in groups at t2 is 



Ni Ni rii 

M'grpitl) = Yl fi'^t = E E I^Ua- 
i i k 



(C2) 
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However, the stellar mass of any particle may differ at t2, so that Mgj.p{t2) 7^ M|^p(ti). The 
difference is, of course, related to the net gain or loss due to star formation and stellar evolution 
mass loss: 

- Mf = M|,p(t2) - M|,p(ii). (C3) 
Similarly, the change in gas mass implies that 

Mf - M'g = M|,^(t2) - M|,^(ii). (C4) 

Since there is no stellar evolution mass loss, M| = 0. The treatment for accreted matter is 
analogous. 
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